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ABSTRACT 


Consider  the  aodel 

x(t)  »  As(t)  +  n(t),  t  *  1 . N 

/v  /v 

where  at  tine  t,  x(t)  is  a  p-coaplex  vector  of  observations,  s(t)  is  a 
q-coaplex  vector  of  unobservable  signals  eaitted  froa  q  sources  and  n(t)  is 
a  noise  p-coaplex  vector  variable.  It  is  assuaed  that  {s(t>}  and  (n(t)} 
are  independent  sequences  with 

E[n(t)n(t)*]  *  <7*1,  E[s(t)s(t  )]  *  T 

/v  /v  a*  a# 

=>  E[x(t)x(t)*]  »  ATA*  +  <r*I 


The  aatrix  A  has  a  special  structure  with  its  k-th  column  ak  of  the  fora 
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where  j  =  V-l,  ^  =  c_1Asin^,  ^  being  the  direction  of  arrival  of  the 

signal  from  the  k-th  source.  The  problem  is  the  estimation  of  q  (the 
number  of  sources)  and  0.’s  (the  directions  of  arrival). 

The  paper  presents  an  information  theoretic  criterion  to  decide  on  q, 
the  number  of  sources,  and  a  spectral  analytic  method  to  estimate  the  r.’s. 

The  proposed  method  is  straight  forward  and  has  some  advantages  over  the 
search  methods  proposed  in  the  literature  for  the  latter  problem.  For 
instance,  our  method  works  even  when  the  signals  are  coherent,  i.e.,  the 
covariance  matrix  T  of  signal  is  singular,  whereas  the  MUSIC  algorithm 
generally  used  in  searching  for  the  estimates  of  is  not  applicable. 

Key  Words  and  Phrases:  Akaike  information  criterion.  Direction  of 
arrival.  General  information  criterion,  Signal  processing. 

AMS  Classification  Index:  Primary  62F12:  Secondary  62H12 


1.  INTRODUCTION 


The  problems  of  estimation  of  the  number  of  signals  and  the  direct 
of  their  arrival  have  been  considered  by  a  number  of  authors  under  the 
model 


x(t)  *  As(t)  +  n(t),  t  *  1,...,N 

M  (V 


(1.1) 


ions 


where  at  time  t, 

x(t)  is  a  p-complex  vector  of  observations  received  by  j  ‘ 

**  I  1 

p-sensors  uniformly  spaced,  j 

s(t)  is  a  q-complex  vector  of  unobservable  signals,  and  j  ; 

j  r 


,D  \v;. 


* 
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n(t)  is  an  additive  p-complex  noise  vector. 

It  is  assuaed  that  {n(t)},  {s(t)},  t  =  1,...,N,  are  independent  sequences 
of  iid  random  variables  such  that 


E(n(t) ]  =  0,  E[n(t)n(t)*]  ■  <r*I 

«v  <v 

E[s( t) ]  =  0,  E[s(t)s(t)*]  =  T  (1.2) 


where  the  star  sign  indicates  the  coaplex  conjugate.  Then  we  have 


E  »  E[x(t)x(t)*]  =  ATA*  +  a1 I. 


(1.3) 


Ve  generally  assuae  that  the  signals  are  not  coherent,  i.e.,  the  matrix  T 
is  positive  definite.  However,  this  condition  is  not  necessary  for  the  new 
method  of  estiaation  we  develop  in  Soction  4  of  this  paper. 


The  pxq  aatrix  A  has  the  special  structure 


As  (a  : . . . :a  ) 

~l  q 

-j  unr,  -jwft(p-i)rv 

a^  =  a(r.)  s  (l,e  0  k,  ...,  e  k) 

Ml  ~  K 


(1.4) 


where  j  s  /If,  s  C-1A  sind^,  c  »  speed  of  propogation,  0^  is  the 

direction  of  arrival  (DOA)  froa  the  k-th  source  and  A  is  the  inter 
sensor  distance.  Ve  nay  take  wQ  to  be  unity  without  loss  of  generality. 

There  are  two  iaportant  probleas  connected  with  the  aodel  (1.1,  1.2). 
One  is  the  estiaation  of  q,  the  number  of  signals  (i.e.,  the  cho'ce  of  the 
aodel),  and  another  is  the  estimation  of  providing  the 
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estimates  of  d , .... ,0  ,  the  directions  of  arrival.  Ve  discuss  these  two 
l  <1 

problems  with  references  to  earlier  work  and  propose  new  solutions. 


2.  INFORMATION  THEORETIC  CRITERIA  FOR  MODEL  SELECTION 

In  view  of  the  assumptions  (1.2),  if  we  assume  that  n(t)  has  a 
p-variate  complex  normal  distribution,  then  the  statistic 

S  =  N-1£  x(t)x(t)*  (2.1) 

has  a  complex  Vishart  distribution  with  N  degrees  of  freedom  and  covariance 

matrix  £  *  aTA*  +  a7 I  as  in  (1.3).  If  the  covariance  matrix  of  the 
signals  is  nonsingular,  then  the  number  of  signals  is  q  *  p-s  where  s 
is  the  multiplicity  of  the  smallest  root  of  £.  Hence  the  problem  of 
estimating  q  can  be  studied  within  the  framework  of  testing  the  equali  ty 
of  a  given  number  of  smallest  eigen  values  of  £,  for  which  a  statisfactory 
solution  exists  (see  for  instance  Anderson  (1963),  Ligget  (1973),  Rao 
(1973,  1983),  Schmidt  (1981),  Tuft  and  Kumaresan  (1980),  Wax,  Shan  and 
Kailath  (1984)). 

The  log  likelihood  based  on  observed  data,  apart  from  an  additive 
constant,  is 


=  -  fUog  |E|  +  tr  FlS]  (2.2) 

where  q  indicates  the  number  of  signals.  The  maximum  of  (2.3)  for  given 
q  is 

^NQ)  "  -  f[j  108  ^i  +  (P“<l)l0S  Vlplq‘  P 


(2.3) 
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where  ^  >  ...  >  are  the  eigen  values  of  S  assumed  to  be  distinct 

(with  probability  1).  Hie  likelihood  ratio  statistic  for  testing  the 
equality  of  the  last  p-q  eigen  values  of  E  is 


%  ¥  +X 

G^q>  »  2(^p,-^q>)  =  N[ (p-q) log  P  -  log($ 


q+1 


••  V1 

(2.4) 


•  a 

which  is  asymptotically  distributed  as  x  on  (p-q)  -  1  degrees  of 

freedom  (d.f.).  Let 

q  »  max  {q:  G^q>  <  Ca) 


m  • 

where  Ca  is  the  upper  tail  crib  point  of  x  on  (p-q)  -  1  d.f.  Then  an 

upper  (l-a)%  confidence  limit  to  q  is  q. 

However,  to  get  a  point  estimate  of  q,  we  can  use  a  general 
information  criterion  (GIC)  developed  by  Zhao,  Krishna! ah  and  Bai  (1986) 

GIC(k)  =  2^k)  -  i/(k)CH  (2.5) 

where  i/(k)  is  the  number  of  free  parameters  in  the  model  and  is  a 

function  of  N  such  that  as  N  -»  « 

iF- 0  ,nd  <2-8> 

In  the  present  problem  i/(k)  *  q(2p-q)  +  1,  and  the  decision  (estimation) 

A 

rule  is  the  choose  q  =  q  where 
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GIC(q)  =  max  GIC(k) .  (2.7) 

k 

It  is  shown  by  Zhao,  Krishnaiah  and  Bai  (1986)  that  under  the  condition 
(2.6),  the  GIC  rule  (2.7)  is  strongly  consistent,  i.e., 

q  -»  q  a.s.  as  n  -»  a.  (2.8) 

Criteria  of  the  kind  (2.7)  for  model  selection  have  been  used  earlier 
with  special  choices: 

Cjj  =  2,  Aka  ike 

a  log  N,  Kashyap,  Schwartz,  Rissanen 
s  c  loglog  n,  for  some  e  >  2,  Hannan  and  Quinn. 

All  these  choices  except  Akaike’s  give  a  consistent  estimate  of  q.  There 
are,  however,  no  adequate  studies  on  the  choice  of  in  finite  samples. 

Rao  and  Vu  (1989)  indicated  an  adaptive  method  for  choosing  in  a 

similar  context.  Further  research  in  this  direction  is  necessary. 


3.  ESTIMATION  OF  DO  A:  NOISE  SPACE  METHOD 


For  a  given  q,  the  eigen  values  of  E  are  of  the  form 

A,  >  ...  >  A  >  A  «...  3  A  (3  <r*).  (3.1) 

1  -  q  q+i  P 

Let  ®!t • • • »Sq.Sq+i* • • ke  the  corresponding  eigen  vectors  and  define 
the  matrices 


(e 


* 


E. 


(£q+l ’ 


.O  . 

~n 


(3.2) 
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The  spaces  spanned  by  the  vectors  in  Eg  and  En  may  be  called  the  signal 
and  noise  eigen  spaces,  respectively.  It  is  easily  seen  that 

a(r)'E„  =  0  <*=>  D(r)  4  a(r)'E Va(r)  *  0  (3.3) 

~  n  ~  n  n~ 

for  r  =  fj,  i  =  l,...,q,  as  defined  in  (1.4).  In  practice  we  have  only 

the  estimate  S  (as  defined  in  (2.1))  of  E.  The  eigen  values  and  eigen 
vectors  of  S 


X 


1 


> 


A 


A 


(3.4) 


provide  consistent  estimates  of  the  corresponding  eigen  values  (3.1)  and 
eigen  vectors  (3.2)  of  the  true  E.  An  estimate  of  En  spanning  the  noise 

eigen  space  is 


n 


(Vi* 


The  function 


a,  ,  ,  . • A  A 

D(r)  x  a(r)  E  E  a(r) 
n  n 


(3.5) 


(3.6) 


which  is  an  estimate  of  D(r)  defined  in  (3.3)  may  not  vanish  for  any  r, 
but  is  likely  to  be  small  in  the  neighborhood  of  fj . rq.  Then  they  may 

be  estimated  by  plotting  1/D(r)  and  locating  the  peaks,  which  is  the 
basis  of  the  MUSIC  algorithm  (see  papers  by  Bienvenu  (1979),  Schmidt  (1981) 
and  Wax,  Shan  and  Kailath  (1984)).  As  the  solutions  are  obtained 
essentially  by  a  search  method,  their  analytical  characterization  is 
complicated  so  that  the  statistical  properties  of  the  estimates  are 
difficult  to  study.  Ve  propose  an  alternative  simple  computational 
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technique  which  provides  definitive  estimates  and  which  enable  us  to  study 
their  asymptotic  properties. 

Ve  note  that  there  exists  a  matrix  G  of  order  pxr  (with  r  written 
for  p-q) 


G  = 


(Gi 


••Sr* 


(3.7) 


where 


G^  is  of  the  form 

£k  =  (0 . 0,gl’ ' • • ’*q+i«0’ 


k  =  1, 


»r, 


(3.8) 


with  the  first  k  and  the  last  (p-k-q-1)  components  as  zeroes 

gq+l  >  °*  Ml  +  ••  +  *q+l*q+l  =  1  *UCh  that  S(r)#G  =  °*  for 
t  =  tx . r^,  and  exp( jrL) , . . . ,exp( jr^)  are  the  solutions  of  the 

polynomial  equation 


Viz<1  +  •••  +  *1  ■  °*  (39) 

Then  the  strategy  is  to  estimate  g  =  (gj , . . . ,gq+1 )  and  use  them  in  (3.9) 
to  obtain  the  roots  which  provide  the  estimates  of 


Since  the  columns  of  G  generate  (i.e.,  provide  a  basis  of)  the  noise 
eigen  space  EQ  of  E,  the  statistical  problem  of  estimating  G  may  be 

thought  of  as  fitting  a  basis  of  the  type  Gj,...,Gr  to  the  estimated 

noise  eigen  space  Er  given  in  (3.5).  The  mathematical  problem  may  be 
formulated  as  that  of  minimizing  the  Euclidean  norm 

l|En-GB|| 


(3.10) 
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with  respect  to  a  matrix  G  of  the  type  (3.7)  and  an  arbitrary  matrix  B  of 
order  rxr. 

As  the  optimization  problem  (3.10)  is  complicated,  we  suggest  some 
alternative  procedures.  By  Householder  transformation  (for  which  reliable 
software  exists),  i.e.,  multiplying  by  an  unitary  matrix  0  of  order  r,  we 

A 

ean  convert  EM  into  the  form 
n 


®n°  -  "  =  <V.-V . V 


(3.11) 


where  u^.  =  (u1>q+i . Sq+i,q+i’° . 0)  ’  *  =  1 . r’  With 

u  .  *  0  (with  probability  1).  le  may  choose  u  ,  as  an  estimate 

q+i ,q+i  '  ~q+l 

of  (gt . gq+1 ) ' ,  solve  the  equation 


B(z)  =  VlQ+l2  +  ••  +  UI  ,q+l  =  0 


(3.12) 


obtain  the  roots  in  the  form 


A 

if. 


A  »'  V 

p^e  ,  k  *  l,...,q 


(3.13) 


A  A 

and  choose  r1>---.*'q  *»  estimates  of  r1...*,fq- 

To  obtain  the  asymptotic  distribution  of  the  estimates,  we  introduce 
some  conditions. 


(At) 

The 

second  moment 

of 

x(t) 

exists . 

(A,) 

The 

fourth  moment 

of 

x(t) 

exists . 
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(Ag)  cov[Re  s ( t ) ]  =  cov[ Im  s(t)]  =  2_1Re  T 

E[(Re  s(t))(Im  s(t))  ]  =  -E[(In  s(t))(Re  s(t))  ]  =  2-1  Im  T 

(A4)  cov[Re  f(t)]  =  cov[Ia  e(t)]  =  2-l<72Ip 
cov[Re  e(t),  Im  f(t)]  =  0 

(A$)  E(Re  fk(t))4  =  E(Im  ek(t))4  =  |<t3,  k=l . p 

E(Re  ck(t))J(Re  ffa(t))3  =  E(Im  fk(t))2(Im  ch(t))2  =  4~V 

k  *  h  =  1 . p 

E(Re  ek(t))2(Im  *h(t))2  =  4“V,  k,h  =  1 . p 

and  all  other  kinds  of  fourth  moments  are  zero. 

The  following  theorems  concerning  the  asymptotic  properties  of  the 
estimators  are  established  in  Bai,  Miao  and  Rao  (1990). 

Theorem  1.  Under  the  condition  (Aj), 


A  $ 


t  =  (r1 . r  )'  t  =  (t1 . r  ) 


(3.14) 


Theorem  2.  Under  the  conditions  (Aj)-(A5)  the  limiting  distribution 


of 


^  +  jTn  =  Mp1-p1 . ?q-pq)  +  jVfT(r1-rl . rQ-rq) 


q  q 


(3.15) 


is  q-variate  complex  normal  with  zero  mean  and  covariance  matrix 


g_1  [<r4r 1  ( A*A)~1r 1  +  <r2rl](G"1)* 


(3.16) 
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whtre 


G  =  diag(D(rl), . . . ,D(rq)) 

D(r,)  =  gfl..  9  (eJ  k  -  eJ  1)eJ  k,  k  ■  1 . q. 

*  q+1  i=l ,  i*k 


The  asymptotic  covariance  matrix  of 


•vv 


alone  is 


2~1Re[G”1  («r4r 1  (A*Ar1r1+<r3r-1  )G_1  ] 


(3.17) 


Remark  1.  Note  that  u  .  . . u  is  a  basis  of  the  eigen  space 

MJ+S+1  ~p 

A  A 

spanned  by  «q+8+1 . ep,  so  that  if  the  number  of  signals  is  q+s,  then 

the  corresponding  g-vector  with  q+s+1  components  is  estimated  by  2q+g+1 • 

In  practice,  it  is  advisable  to  estimate  the  r  parameters  starting  with  a 
somewhat  higher  value  than  what  is  indicated  by  the  GIC  criterion,  study 
the  configuration  of  the  roots  as  done  by  Tuft  and  Eumaresan  (1980)  and 
take  a  final  decision  on  the  choice  of  q  and  estimate  the  r  parameters. 
The  computational  algorithm  outlined  above  is  ideally  suited  to  such  an 
analysis . 

Remark  2.  We  can  by  postmultiplying  H  in  (3.11)  by  a  matrix  B 
reduce  it  to  the  form 


HB 


(u  V 
vJ5q+l*Zq+2» 


(3.18) 


where  =  (0, . . .  ,0,v1  q+J , . . .  ,y  j  ,0, . . .  ,0)  with  i-1  zeroes  in  the 

beginning  and  ||Vq+jl!  =  1.  Then  a  possibly  improved  estimate  of  gj  in 

(*i . «,«•  »• 
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A 


Uj.q+l 


+ 


vj,q+a 


+ 


Vj  p),  j  =  1 . q+1.  (3.19) 


The  actual  improvement  in  the  estimate  (3.19)  needs  to  be  studied. 


Remark  3.  There  is  a  unique  vector  of  the  type  u^+1  in  the  space 

generated  by  the  columns  of  E^.  This  can  be  obtained  in  a  simple  way, 
without  going  through  Householder  transformation,  as  shown  below.  We 

A 

partition  in  the  form 


En 

px(p-q) 


r\ 

%1 

(q+i)x(p-q) 

A 

h, 

Up-q-l)x(p-q)J 


and  solve  the  equation 


,b,) 

p-q 


(3.20) 


(3.21) 


The  solution  of  (3.21)  is  unique  (with  probability  1)  apart  from  a 
scalar  multiplier.  Then 


t  *  <*! 


apart  from  a  scalar  multiplier. 


4.  ESTIMATION  OF  DO  A:  DIRECT  APPROACH 


Let  us  write  G. 


as  defined  in  (3.8)  in  the  partitioned  form 
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G-  =  (  0  :  g  :  0  ) 

lX(i-l)  lX(q+l )  lX(p-q-l) 


(4.1) 


where 


«-<*, . *<i«>0'  U*H  •  i.  ‘i  =  «- 


(4.2) 


Since  Gj  is  an  eigen  vector  of  E  *  ATA*  +  a7l  corresponding  to  the 

smallest  eigen  value  <r3  for  each  i,  it  follows  that  for  any  chosen 
non-negative  , . . .  ,1/^^ 


pMJ  ^  g  P- (J  A  <f| 

z*'kSk£Sk-*<2  •'Ao,«  4  *  V 


(4.3) 


attains  the  ainimum  value  ai{vl  +  ...  +  at  S  satisfying  (4.2), 

where  £jeQ  is  the  matrix  obtained  from  E  by  retaining  only  the  k-th  to 
(k+q)_th  columns  and  rows,  then  g  may  be  estimated  by  minimizing 


>-q  *  . 

E  Sk  "  Sk’  subiect  *°  (4*2). 


(4.4) 


Writing  for  the  (q+l)x(q+l)  matrix  obtained  by  retaining  only  the 

k-th  to  (k+q)-th  columns  and  rows  in  £,  the  expression  (4.4)  can  be 


written  as 


**(T,Ao)*  * 


(4.5) 
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so  that  g,  the  estimate  of  g,  is  the  eigen  vector  of  corresponding 

to  its  smallest  eigen  value.  Then  we  solve  the  polynomial  equation 


A  k 
*q+lZ  + 


A 

+  s. 


(4.6) 


as  in  Section  3.  Writing  the  roots  in  the  form 

A  A 

>  0,  e  (0,2*),  we  obtain  the  estimate  of 


=  Pk  exp(jrk), 

A 

rk  as  rk’  k  =  1, . . .  ,q. 


Note  that  in  (4.3) 


p— q  p— q  i  ,  p— q 

*  E  ^k^ko  ~  ^  1  +  (<rz  E  vk)Iq+1  (4.7) 

so  that  even  if  T  is  singular,  when  p  is  large  enough,  the  first  term 
in  (4.7)  can  be  of  rank  q.  The  classical  MUSIC  algorithm  fails  when  T  is 
singular,  but  the  above  method  works. 


A  A 

In  order  to  study  the  limiting  distribution  of  p^,  r k,  k  *  1 . q, 

we  make  the  following  notations: 


jrk  i<Tv 

Sko  =  (1,e  . 8  (q+l)-vector. 

Ao  *  (Sio)'--"Sqo)*  (<*+1)x<*  Mtrix. 
jr  jr 

B  *  diag(e  . e  H),  qxq  diagonal  matrix. 

Theorem  1.  Under  the  condition  (A}) 

(4.8) 

(4.9) 


l  *  (r! . V 


q' 

A  A  # 

(P| » • • • ) 


a. s. 


I  *  . V 

£  *  (1 . 1)'. 


t» 
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Theorem  2.  Under  the  condition  (Aj),  the  joint  distribution  of 

A 

+  jT„,  where 
n 

K  =  ^(jri . aq- 1) 

A  /—A  A 

Tn  *  vN(rl-rl . r  -r  ), 

tends  to  that  of 


l 


kBk"lr  B  k“1)“1(A*A0)“1A*(Piqi/k(Rjk)+A0Bk~l 


Rjk>)) 


(4.10) 


where  (Rjk,,R‘k)) 


has  a  Multivariate  coaplex  noraal  distribution  and 


q+l  jr  (k-l) 

Q  =  diag(  E  (k-l)g,e  1 
fc=l  * 


q+i 

,  E  (k-l)g.e  q 
k»i  k 


(k-i  > 


). 


(4.11) 


Let  R  be  a  pxp  Gaussian  random  aatrix  with  the  underlying  variance 

<r4  and  V  be  a  qxp  randoa  aatrix  whose  eoluans  are  iid  aultivariate 

* 

complex  noraal  variables  with  aean  zero  and  covariance  aatrix  a  V. 


Theorea  3.  Under  the  conditions  (A3)-(Aj),  the  randoa  variables  R‘k 
and  R*k)  in  (4.10)  can  be  characterized  as  follows: 

(i)  Rjk>  is  the  subaatrix  of  V  foraed  by  its  (k+l)-th  to 
(k+q+l)-th  columns. 


(ii)  Rj  1  is  the  submatrix  formed  by  the  (k+l)-th  to 
(k+q+l)-th  comumns  and  rows  of  R. 


(iii)  (Rj* 1 , . . . ,RjP“<l) )  and  (R^1 ’ , . . . ,Rjp-q) )  are  independently 


<p-<n 


distributed. 


From  Theorems  2  and  3,  it  follows  that  the  limiting  distribution  of 


A  A 


ip  +  jTQ  is  q-variate  complex  normal  with  mean  zero  and  covariance  matrix 


Q’1(Piqi/lrBk“1r  b  k"1)“lv(p?i/vBk_1r  b  k-1)Q-1 

l  *  k*i  * 


(4.12) 


V  =  (aXj-'aJe  E  ‘V'sFks*ks<7VVo)"1  +  E  E  gk8  (4.13) 


gks  58  E 


Ka  »  *  1  »  •  • 

-k  *1  "a 


(4.14) 


if  k  >  s 

Pfcs  *  l(J8”V  if  k  <  s 


(4.15) 


where  J  is  a  (p-q)x(p-q)  matrix  with  only  the  entries  in  the  diagonal 
above  the  main  diagonal  as  unities  and  rest  as  zeroes.  The  asymptotic 

A  _i 

covariance  of  T  alone  is  2  times  the  real  part  of  (4.12) 

Remark  4.1.  It  is  seen  that  the  expressions  for  the  covariance 

A  A 

matrices  of  the  estimates  r  and  g  given  in  (3.16)  and  (4.12)  are  very 

complicated.  Perhaps,  here  we  have  a  case  for  the  bootstrap  method  for 
estimating  the  covariance  matrix.  The  numerical  computations  involved  in 
the  bootstrap  method  will  be  considered  in  a  later  communication. 
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Remark  4.2.  It  is  seen  from  (1.1)  that  the  conjugate  vector  x(t) 
has  the  aodel 


x(t)  »  A  s(t)  +  n(t),  t  =  1 . N. 


(4.16) 


with  the  covariance  matrix  of  x(t)  as 


S  =  A  T  A*  +  <r2I. 


(4.17) 


Then 


Si  *  (0 . 0,gq+i'  *  ‘  *  ,gl*° . 0)  (4.18) 

which  is  siailar  to  6.  of  (3.8)  with  the  components  of  the  same  g  vector 
written  in  the  reverse  order,  is  an  eigen  vector  corresponding  to  the 

smallest  eigen  value  of  2.  Ve  can  estimate  gq+x**--»gi  using  the  model 

(4.16)  in  the  same  way  as  gx>'**«gq+1  **  estimated  using  the  model  (1.1). 

Ve  have  a  developed  method  of  combining  the  information  proved  by  the 
models  (1.1)  and  (4.16)  in  obtaining  an  efficient  estimate  of  glt...,gq+1, 

which  will  be  given  in  a  later  communication. 
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